#define kirby2_N       151
#define kirby2_P       5

static double kirby2_x0a[kirby2_P] = { 2.0, -0.1, 0.003, -0.001, 0.00001 };
static double kirby2_x0b[kirby2_P] = { 1.5, -0.15, 0.0025, -0.0015, 0.00002 }; 
static double kirby2_x[kirby2_P] = { 1.6745063063E+00, -1.3927397867E-01, 2.5961181191E-03,
                                    -1.7241811870E-03,  2.1664802578E-05 }; 

static double kirby2_epsrel = 1.0e-5;

static double kirby2_sigma[kirby2_P] = {
  8.7989634338E-02, 4.1182041386E-03, 4.1856520458E-05,
  5.8931897355E-05, 2.0129761919E-07
};

static double kirby2_F1[kirby2_N] = {
       0.0082E0, 0.0112E0, 0.0149E0, 0.0198E0, 0.0248E0,
       0.0324E0, 0.0420E0, 0.0549E0, 0.0719E0, 0.0963E0,
       0.1291E0, 0.1710E0, 0.2314E0, 0.3227E0, 0.4809E0,
       0.7084E0, 1.0220E0, 1.4580E0, 1.9520E0, 2.5410E0,
       3.2230E0, 3.9990E0, 4.8520E0, 5.7320E0, 6.7270E0,
       7.8350E0, 9.0250E0, 10.2670E0, 11.5780E0, 12.9440E0,
      14.3770E0, 15.8560E0, 17.3310E0, 18.8850E0, 20.5750E0,
      22.3200E0, 22.3030E0, 23.4600E0, 24.0600E0, 25.2720E0,
      25.8530E0, 27.1100E0, 27.6580E0, 28.9240E0, 29.5110E0,
      30.7100E0, 31.3500E0, 32.5200E0, 33.2300E0, 34.3300E0,
      35.0600E0, 36.1700E0, 36.8400E0, 38.0100E0, 38.6700E0,
      39.8700E0, 40.0300E0, 40.5000E0, 41.3700E0, 41.6700E0,
      42.3100E0, 42.7300E0, 43.4600E0, 44.1400E0, 44.5500E0,
      45.2200E0, 45.9200E0, 46.3000E0, 47.0000E0, 47.6800E0,
      48.0600E0, 48.7400E0, 49.4100E0, 49.7600E0, 50.4300E0,
      51.1100E0, 51.5000E0, 52.1200E0, 52.7600E0, 53.1800E0,
      53.7800E0, 54.4600E0, 54.8300E0, 55.4000E0, 56.4300E0,
      57.0300E0, 58.0000E0, 58.6100E0, 59.5800E0, 60.1100E0,
      61.1000E0, 61.6500E0, 62.5900E0, 63.1200E0, 64.0300E0,
      64.6200E0, 65.4900E0, 66.0300E0, 66.8900E0, 67.4200E0,
      68.2300E0, 68.7700E0, 69.5900E0, 70.1100E0, 70.8600E0,
      71.4300E0, 72.1600E0, 72.7000E0, 73.4000E0, 73.9300E0,
      74.6000E0, 75.1600E0, 75.8200E0, 76.3400E0, 76.9800E0,
      77.4800E0, 78.0800E0, 78.6000E0, 79.1700E0, 79.6200E0,
      79.8800E0, 80.1900E0, 80.6600E0, 81.2200E0, 81.6600E0,
      82.1600E0, 82.5900E0, 83.1400E0, 83.5000E0, 84.0000E0,
      84.4000E0, 84.8900E0, 85.2600E0, 85.7400E0, 86.0700E0,
      86.5400E0, 86.8900E0, 87.3200E0, 87.6500E0, 88.1000E0,
      88.4300E0, 88.8300E0, 89.1200E0, 89.5400E0, 89.8500E0,
      90.2500E0, 90.5500E0, 90.9300E0, 91.2000E0, 91.5500E0,
      92.2000E0
};


static double kirby2_F0[kirby2_N] = {
    9.65E0, 10.74E0, 11.81E0, 12.88E0, 14.06E0,
   15.28E0, 16.63E0, 18.19E0, 19.88E0, 21.84E0,
   24.00E0, 26.25E0, 28.86E0, 31.85E0, 35.79E0,
   40.18E0, 44.74E0, 49.53E0, 53.94E0, 58.29E0,
   62.63E0, 67.03E0, 71.25E0, 75.22E0, 79.33E0,
   83.56E0, 87.75E0, 91.93E0, 96.10E0, 100.28E0,
  104.46E0, 108.66E0, 112.71E0, 116.88E0, 121.33E0,
  125.79E0, 125.79E0, 128.74E0, 130.27E0, 133.33E0,
  134.79E0, 137.93E0, 139.33E0, 142.46E0, 143.90E0,
  146.91E0, 148.51E0, 151.41E0, 153.17E0, 155.97E0,
  157.76E0, 160.56E0, 162.30E0, 165.21E0, 166.90E0,
  169.92E0, 170.32E0, 171.54E0, 173.79E0, 174.57E0,
  176.25E0, 177.34E0, 179.19E0, 181.02E0, 182.08E0,
  183.88E0, 185.75E0, 186.80E0, 188.63E0, 190.45E0,
  191.48E0, 193.35E0, 195.22E0, 196.23E0, 198.05E0,
  199.97E0, 201.06E0, 202.83E0, 204.69E0, 205.86E0,
  207.58E0, 209.50E0, 210.65E0, 212.33E0, 215.43E0,
  217.16E0, 220.21E0, 221.98E0, 225.06E0, 226.79E0,
  229.92E0, 231.69E0, 234.77E0, 236.60E0, 239.63E0,
  241.50E0, 244.48E0, 246.40E0, 249.35E0, 251.32E0,
  254.22E0, 256.24E0, 259.11E0, 261.18E0, 264.02E0,
  266.13E0, 268.94E0, 271.09E0, 273.87E0, 276.08E0,
  278.83E0, 281.08E0, 283.81E0, 286.11E0, 288.81E0,
  291.08E0, 293.75E0, 295.99E0, 298.64E0, 300.84E0,
  302.02E0, 303.48E0, 305.65E0, 308.27E0, 310.41E0,
  313.01E0, 315.12E0, 317.71E0, 319.79E0, 322.36E0,
  324.42E0, 326.98E0, 329.01E0, 331.56E0, 333.56E0,
  336.10E0, 338.08E0, 340.60E0, 342.57E0, 345.08E0,
  347.02E0, 349.52E0, 351.44E0, 353.93E0, 355.83E0,
  358.32E0, 360.20E0, 362.67E0, 364.53E0, 367.00E0,
  371.30E0
};

static void
kirby2_checksol(const double x[], const double sumsq,
                const double epsrel, const char *sname,
                const char *pname)
{
  size_t i;
  const double sumsq_exact = 3.9050739624E+00;

  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
               sname, pname);

  for (i = 0; i < kirby2_P; ++i)
    {
      gsl_test_rel(x[i], kirby2_x[i], epsrel, "%s/%s i=%zu",
                   sname, pname, i);
    }
}

static int
kirby2_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double b[kirby2_P];
  size_t i;

  for (i = 0; i < kirby2_P; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < kirby2_N; i++)
    {
      double t = kirby2_F0[i];
      double y = ((b[0] + t* (b[1]  + t * b[2]))
                  / (1 + t*(b[3]  + t *b[4])));
      gsl_vector_set (f, i, kirby2_F1[i] - y);
    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;
}

static int
kirby2_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double b[kirby2_P];
  size_t i;

  for (i = 0; i < kirby2_P; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < kirby2_N; i++)
    {
      double t = kirby2_F0[i];
      double u = (b[0] + t*(b[1] + t*b[2]));
      double v = (1 + t*(b[3] + t*b[4]));
      gsl_matrix_set (df, i, 0, -1/v);
      gsl_matrix_set (df, i, 1, -t/v);
      gsl_matrix_set (df, i, 2, -t*t/v);
      gsl_matrix_set (df, i, 3, t*u/(v*v));
      gsl_matrix_set (df, i, 4, t*t*u/(v*v));
    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;
}

static int
kirby2_fvv (const gsl_vector * x, const gsl_vector * v,
            void *params, gsl_vector * fvv)
{
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
  double x4 = gsl_vector_get(x, 3);
  double x5 = gsl_vector_get(x, 4);
  double v1 = gsl_vector_get(v, 0);
  double v2 = gsl_vector_get(v, 1);
  double v3 = gsl_vector_get(v, 2);
  double v4 = gsl_vector_get(v, 3);
  double v5 = gsl_vector_get(v, 4);
  size_t i;

  for (i = 0; i < kirby2_N; i++)
    {
      double t = kirby2_F0[i];
      double term1 = 1.0 + t*(x4 + t*x5);

      gsl_vector_set(fvv, i, -2*t/pow(term1, 3.0) *
                             (v4 + t*v5) *
                             (-t*(-v2 + v4*x1 +
                                  t*(-v3 + v5*x1 + v4*x2 + t*v5*x2 + t*(v4+t*v5)*x3)) +
                             t*t*(v2 + t*v3)*(x4 + t*x5) + v1*(1.0 + t*(x4 + t*x5))));
    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;
}

static gsl_multifit_nlinear_fdf kirby2_func =
{
  kirby2_f,
  kirby2_df,
  kirby2_fvv,
  kirby2_N,
  kirby2_P,
  NULL,
  0,
  0,
  0
};

static test_fdf_problem kirby2a_problem =
{
  "nist-kirby2a",
  kirby2_x0a,
  NULL,
  kirby2_sigma,
  &kirby2_epsrel,
  &kirby2_checksol,
  &kirby2_func
};

static test_fdf_problem kirby2b_problem =
{
  "nist-kirby2b",
  kirby2_x0b,
  NULL,
  kirby2_sigma,
  &kirby2_epsrel,
  &kirby2_checksol,
  &kirby2_func
};
